Spotify Data Analysis and Descriptive Statistics

1. Data Set Description For this project, our team decided to take a data-driven approach to evaluate music. We are using the Spotify Dataset from 1921 to 2020 that consists of over 160,000 different tracks. In order to stay consistent in the evaluation, all data was sourced from the Spotify Web API. The following set of data is combination of Primary, Numerical, Dummy(binary), and Categorical data. This allows us to explore different types of models and draw different insights. Here are all the variables that is included in the data set:

Primary - id: this an unique key comprised of numbers and characters that is assigned to each track generated by Spotify

Numerical - acousticness: range from 0(LOW) to 1(HIGH) - danceability: range from 0(LOW) to 1(HIGH) - energy: range from 0(LOW) to 1(HIGH) - duration_ms: majority range from 200,000 to 300,000 - instrumentalness: range from 0(LOW) to 1(HIGH) - valence: range from 0(LOW) to 1(HIGH) - popularity: range from 0(LOW) to 100(HIGH)* - tempo: majority range from 50(LOW) to 150(high) - liveness: range from 0(LOW) to 1(HIGH) - loudness majority range from -60 to 0 - speechiness: range from 0(LOW) to 1(HIGH) - year: range from 1921 to 2020

Dummy - mode: 0 represents minor and 1 represents major - explicit: 0 represents no explicit content and 1 represents explicit content

Categorical - key: this consists of all different music keys on onctave encoded from 0 to 11 (i.e. C = 0, C# = 1, etc…) - artists: the artist of the track - release_date: the date of release in yyyy-mm-dd format - name: the name of the track

*With our approach on evaluating different tracks, we decided to use popularity as our main dependant variable.

options(repr.matrix.max.rows=100, repr.matrix.max.cols=20)

Import all Libraries

# Libraries
options(warn=-1)
library(ggplot2)
library(dplyr)
library(plotly)
library(hrbrthemes)
library(forecast)
library(xts)
library(Metrics)
library(psych)
library(dygraphs)
library(GGally)
library(tidyverse)
library(tidyquant)  
library(cranlogs)   
library(corrr)      
library(cowplot)  
library(Metrics)

Annual data for each feature

df <- read.csv("./data/data.csv")
q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness")
for(i in 1921:2020){
  annual_val = subset(df, year==i)

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)


#acoustics = ts(feature_quantiles[,,'acousticness'], start=1921, end=2020, frequency=1)
#ggplot()+geom_line(data=acoustics[,"Median"], mapping=aes(x=1921:2020, y=as.numeric(acoustics[,"Median"]))) +  geom_point() + geom_errorbar(aes(x = 1921:2020,ymin=as.numeric(acoustics[,"Q1"]), ymax=as.numeric(acoustics[,"Q3"])), )

#geom_line(mapping=aes(x=1921:2020, y = as.numeric(ts(feature_quantiles[,,'liveness'], start=1921, end=2020, frequency=1)[,"Median"])))

#acoustics

plot(1930,1, xlim = c(1921,2020), ylim = range(0,2), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])
  # Normalize it????? / make it a function of its own mean
  
  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5), legend=feature_names)

popular

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  
  
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  annual_val = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    q1s <- rbind(q1s, c(i, name, vals[1]))
    medians <- rbind(medians, c(i, name, vals[2]))
    q3s <- rbind(q3s, c(i, name, vals[3]))
  }
}

colnames(q1s) <-c("Year", "Feature", "Value")
colnames(medians) <-c("Year", "Feature", "Value")
colnames(q3s) <-c("Year", "Feature", "Value")

split_q1s <- split(q1s, q1s$Feature)
split_medians <- split(medians, medians$Feature)
split_q3s <- split(q3s, q3s$Feature)

feature_quantiles = array(0, c(100,3,length(feature_names)))

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,1,i] = split_q1s[[i]][[3]][j]
    feature_quantiles[j,2,i] = split_medians[[i]][[3]][j]
    feature_quantiles[j,3,i] = split_q3s[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    c("Q1", "Median", "Q3"),
                                    feature_names)


#acoustics = ts(feature_quantiles[,,'acousticness'], start=1921, end=2020, frequency=1)
#ggplot()+geom_line(data=acoustics[,"Median"], mapping=aes(x=1921:2020, y=as.numeric(acoustics[,"Median"]))) +  geom_point() + geom_errorbar(aes(x = 1921:2020,ymin=as.numeric(acoustics[,"Q1"]), ymax=as.numeric(acoustics[,"Q3"])), )

#geom_line(mapping=aes(x=1921:2020, y = as.numeric(ts(feature_quantiles[,,'liveness'], start=1921, end=2020, frequency=1)[,"Median"])))

#acoustics

plot(1930,1, xlim = c(1921,2020), ylim = range(0,2), xlab="Year", ylab="Range")
for(i in 1:length(feature_names)){
  name = feature_names[i]
  medians = as.numeric(feature_quantiles[,,name][,"Median"])
  # Normalize it????? / make it a function of its own mean
  
  avg = mean(medians)
  medians = medians / avg
  
  lines(x=1921:2020, y=medians, col=i)
}

legend("topright", lty=1, col=c(1,2,3,4,5), legend=feature_names)

Now trying to plot absolute values

# Get the list of features I need to work with

library(Metrics)

q1s = data.frame(matrix(ncol=3, nrow=0))
medians = data.frame(matrix(ncol=3, nrow=0))
best_medians = data.frame(matrix(ncol=3, nrow=0))
q3s = data.frame(matrix(ncol=3, nrow=0))
feature_names = sort(c("acousticness", "danceability", "instrumentalness", "energy",
                  "duration_ms", "valence", "tempo", "liveness", "loudness", "speechiness"))
for(i in 1921:2020){
  annual_val = subset(df, year==i)
  top_10 = floor(nrow(annual_val)*0.1)
  annual_val = annual_val[order(annual_val$popularity, decreasing=TRUE),]
  best_ones = annual_val[1:top_10,]

  for(name in feature_names){
    vals = quantile(annual_val[,name])[2:4]
    best_vals = quantile(best_ones[,name])[2:4]
    best_medians <- rbind(best_medians, c(i, name, best_vals[2]))
    medians <- rbind(medians, c(i, name, vals[2]))
  }
}

colnames(medians) <-c("Year", "Feature", "Value")
colnames(best_medians) <- c("Year", "Feature", "Value")

split_medians <- split(medians, medians$Feature)
split_best_medians <- split(best_medians, best_medians$Feature)

feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)
best_feature_quantiles = matrix(0, ncol=length(feature_names), nrow=100)

for(i in 1:length(feature_names)){
  for(j in 1:100){
    feature_quantiles[j,i] <- split_medians[[i]][[3]][j]
    best_feature_quantiles[j,i] = split_best_medians[[i]][[3]][j]
  }
}

dimnames(feature_quantiles) <- list(1921:2020,
                                    feature_names)
dimnames(best_feature_quantiles) <- list(1921:2020,
                                    feature_names)




for(i in 1:length(feature_names)){
  name = feature_names[i]
  scores = as.numeric(feature_quantiles[,i])
  best_scores = as.numeric(best_feature_quantiles[,i])
  mape = mape(best_scores, scores)
  mape = round(mape * 100,3)
  r  = cor(scores, best_scores)
  pval = t.test(scores, )
  # Normalize it????? / make it a function of its own mean
  
  smallest = min(min(scores), min(best_scores))
  biggest = max(max(scores), max(best_scores))
  
  plot(x=1921:2020, y=scores, col=1, type="l", ylim=c(smallest, biggest), main=name, sub=paste("MAPE: ", mape, "%;", "R:", r), xlab="Year", ylab="Median Value")
  lines(x=1921:2020, y=best_scores, col=2)
  legend("topright", lty=1, col=c(1,2), legend=c("Overall", "Most Popular"))
  title()
}

NA

Train time series models to forecast the future models. Set up the libraries and the training/testing amounts

library("smooth")
library("forecast")
library("nnfor")
training.percent = 0.95
nTrain = 100*training.percent
nTest = 100*(1-training.percent)

Accousticness Models

acousticness = ts(as.numeric(feature_quantiles[,"acousticness"]) , start=1921)
acousticness.train = subset(acousticness, start=1, end=nTrain)
acousticness.test = subset(acousticness, start = (nTrain+1), end =(nTrain+nTest))

#ses
acousticness.train.ses <- ses(acousticness.train, h=nTest)
acousticness.ses.mape = mape(acousticness.train.ses$mean, acousticness.test)
plot(acousticness.train.ses, main=paste("Acousticness", "SES"), sub=paste("MAPE:", round(acousticness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#SARIMA
acousticness.sarima.model <- auto.arima(acousticness.train)
acousticness.sarima <- forecast(acousticness.sarima.model, h=nTest)
acousticness.sarima.mape = mape(acousticness.sarima$mean, acousticness.test)
plot(acousticness.sarima, main=paste("Acousticness", "SARIMA"), sub=paste("MAPE:", round(acousticness.sarima.mape*100,3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#Neural Network
acousticness.train.nn <- elm(acousticness.train)
acousticness.nn.forecast <- forecast(acousticness.train.nn, h=nTest)
acousticness.nn.mape = mape(acousticness.nn.forecast$mean, acousticness.test)
plot(acousticness, main=paste("Acousticness", "NN"), sub=paste("MAPE:", round(acousticness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(acousticness.nn.forecast$mean, start=1921+nTrain), col=2)

Energy

energy = ts(as.numeric(feature_quantiles[,"energy"]) , start=1921)
energy.train = subset(energy, start=1, end=nTrain)
energy.test = subset(energy, start = (nTrain+1), end =(nTrain+nTest))

#SES
energy.train.ses <- ses(energy.train, h=nTest)
energy.ses.mape = mape(energy.train.ses$mean, energy.test)
plot(energy.train.ses, main=paste("Energy", "SES"), sub=paste("MAPE:", round(energy.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#SARIMA
energy.sarima.model <- auto.arima(energy.train)
energy.sarima <- forecast(energy.sarima.model, h=nTest)
energy.sarima.mape = mape(energy.sarima$mean, energy.test)
plot(energy.sarima, main=paste("Energy", "SARIMA"), sub=paste("MAPE:", round(energy.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#Neural Network
energy.train.nn <- elm(energy.train)
energy.nn.forecast <- forecast(energy.train.nn, h=nTest)
energy.nn.mape = mape(energy.nn.forecast$mean, energy.test)
plot(energy, main=paste("Energy", "NN"), sub=paste("MAPE:", round(energy.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(energy.nn.forecast$mean, start=1921+nTrain), col=2)

Danceability

library("smooth")
library("forecast")
library(nnfor)
training.percent = 0.95
nTrain = 100*training.percent
nTest = 100*(1-training.percent)


#Acousticness
acousticness = ts(as.numeric(feature_quantiles[,"acousticness"]) , start=1921)
acousticness.train = subset(acousticness, start=1, end=nTrain)
acousticness.test = subset(acousticness, start = (nTrain+1), end =(nTrain+nTest))

#ses
acousticness.train.ses <- ses(acousticness.train, h=nTest)
acousticness.ses.mape = mape(acousticness.train.ses$mean, acousticness.test)
plot(acousticness.train.ses, main=paste("Acousticness", "SES"), sub=paste("MAPE:", round(acousticness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#SARIMA
acousticness.sarima.model <- auto.arima(acousticness.train)
acousticness.sarima <- forecast(acousticness.sarima.model, h=nTest)
acousticness.sarima.mape = mape(acousticness.sarima$mean, acousticness.test)
plot(acousticness.sarima, main=paste("Acousticness", "SARIMA"), sub=paste("MAPE:", round(acousticness.sarima.mape*100,3), "%"), xlab="Year", ylab="Median Value")
lines(acousticness)

#Neural Network
acousticness.train.nn <- elm(acousticness.train)
acousticness.nn.forecast <- forecast(acousticness.train.nn, h=nTest)
acousticness.nn.mape = mape(acousticness.nn.forecast$mean, acousticness.test)
plot(acousticness, main=paste("Acousticness", "NN"), sub=paste("MAPE:", round(acousticness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(acousticness.nn.forecast$mean, start=1921+nTrain), col=2)

#Energy
energy = ts(as.numeric(feature_quantiles[,"energy"]) , start=1921)
energy.train = subset(energy, start=1, end=nTrain)
energy.test = subset(energy, start = (nTrain+1), end =(nTrain+nTest))

#SES
energy.train.ses <- ses(energy.train, h=nTest)
energy.ses.mape = mape(energy.train.ses$mean, energy.test)
plot(energy.train.ses, main=paste("Energy", "SES"), sub=paste("MAPE:", round(energy.ses.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#SARIMA
energy.sarima.model <- auto.arima(energy.train)
energy.sarima <- forecast(energy.sarima.model, h=nTest)
energy.sarima.mape = mape(energy.sarima$mean, energy.test)
plot(energy.sarima, main=paste("Energy", "SARIMA"), sub=paste("MAPE:", round(energy.sarima.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(energy)

#Neural Network
energy.train.nn <- elm(energy.train)
energy.nn.forecast <- forecast(energy.train.nn, h=nTest)
energy.nn.mape = mape(energy.nn.forecast$mean, energy.test)
plot(energy, main=paste("Energy", "NN"), sub=paste("MAPE:", round(energy.nn.mape* 100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(energy.nn.forecast$mean, start=1921+nTrain), col=2)

#Duration
duration = ts(as.numeric(feature_quantiles[,"duration_ms"]) , start=1921)
duration.train = subset(duration, start=1, end=nTrain)
duration.test = subset(duration, start = (nTrain+1), end =(nTrain+nTest))

#SES
duration.train.ses <- ses(duration.train, h=nTest)
duration.ses.mape = mape(duration.train.ses$mean, duration.test)
plot(duration.train.ses, main=paste("Duration", "SES"), sub=paste("MAPE:",round(duration.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#SARIMA
duration.sarima.model <- auto.arima(duration.train)
duration.sarima <- forecast(duration.sarima.model, h=nTest)
duration.sarima.mape = mape(duration.sarima$mean, duration.test)
plot(duration.sarima, main=paste("Duration", "SARIMA"), sub=paste("MAPE:", round(duration.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#Neural Network
duration.train.nn <- elm(duration.train)
duration.nn.forecast <- forecast(duration.train.nn, h=nTest)
duration.nn.mape = mape(duration.nn.forecast$mean, duration.test)
plot(duration, main=paste("Duration", "NN"), sub=paste("MAPE:", round(duration.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(duration.nn.forecast$mean, start=1921+nTrain), col=2)

Duration

duration = ts(as.numeric(feature_quantiles[,"duration_ms"]) , start=1921)
duration.train = subset(duration, start=1, end=nTrain)
duration.test = subset(duration, start = (nTrain+1), end =(nTrain+nTest))

#SES
duration.train.ses <- ses(duration.train, h=nTest)
duration.ses.mape = mape(duration.train.ses$mean, duration.test)
plot(duration.train.ses, main=paste("Duration", "SES"), sub=paste("MAPE:",round(duration.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#SARIMA
duration.sarima.model <- auto.arima(duration.train)
duration.sarima <- forecast(duration.sarima.model, h=nTest)
duration.sarima.mape = mape(duration.sarima$mean, duration.test)
plot(duration.sarima, main=paste("Duration", "SARIMA"), sub=paste("MAPE:", round(duration.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(duration)

#Neural Network
duration.train.nn <- elm(duration.train)
duration.nn.forecast <- forecast(duration.train.nn, h=nTest)
duration.nn.mape = mape(duration.nn.forecast$mean, duration.test)
plot(duration, main=paste("Duration", "NN"), sub=paste("MAPE:", round(duration.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(duration.nn.forecast$mean, start=1921+nTrain), col=2)

Valence

valence = ts(as.numeric(feature_quantiles[,"valence"]) , start=1921)
valence.train = subset(valence, start=1, end=nTrain)
valence.test = subset(valence, start = (nTrain+1), end =(nTrain+nTest))

#SES
valence.train.ses <- ses(valence.train, h=nTest)
valence.ses.mape = mape(valence.train.ses$mean, valence.test)
plot(valence.train.ses, main=paste("valence", "SES"), sub=paste("MAPE:",round(valence.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#SARIMA
valence.sarima.model <- auto.arima(valence.train)
valence.sarima <- forecast(valence.sarima.model, h=nTest)
valence.sarima.mape = mape(valence.sarima$mean, valence.test)
plot(valence.sarima, main=paste("valence", "SARIMA"), sub=paste("MAPE:", round(valence.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(valence)

#Neural Network
valence.train.nn <- elm(valence.train)
valence.nn.forecast <- forecast(valence.train.nn, h=nTest)
valence.nn.mape = mape(valence.nn.forecast$mean, valence.test)
plot(valence, main=paste("valence", "NN"), sub=paste("MAPE:", round(valence.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(valence.nn.forecast$mean, start=1921+nTrain), col=2)

Tempo

tempo = ts(as.numeric(feature_quantiles[,"tempo"]) , start=1921)
tempo.train = subset(tempo, start=1, end=nTrain)
tempo.test = subset(tempo, start = (nTrain+1), end =(nTrain+nTest))

#SES
tempo.train.ses <- ses(tempo.train, h=nTest)
tempo.ses.mape = mape(tempo.train.ses$mean, tempo.test)
plot(tempo.train.ses, main=paste("tempo", "SES"), sub=paste("MAPE:",round(tempo.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#SARIMA
tempo.sarima.model <- auto.arima(tempo.train)
tempo.sarima <- forecast(tempo.sarima.model, h=nTest)
tempo.sarima.mape = mape(tempo.sarima$mean, tempo.test)
plot(tempo.sarima, main=paste("tempo", "SARIMA"), sub=paste("MAPE:", round(tempo.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(tempo)

#Neural Network
tempo.train.nn <- elm(tempo.train)
tempo.nn.forecast <- forecast(tempo.train.nn, h=nTest)
tempo.nn.mape = mape(tempo.nn.forecast$mean, tempo.test)
plot(tempo, main=paste("tempo", "NN"), sub=paste("MAPE:", round(tempo.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(tempo.nn.forecast$mean, start=1921+nTrain), col=2)

Liveness

liveness = ts(as.numeric(feature_quantiles[,"liveness"]) , start=1921)
liveness.train = subset(liveness, start=1, end=nTrain)
liveness.test = subset(liveness, start = (nTrain+1), end =(nTrain+nTest))

#SES
liveness.train.ses <- ses(liveness.train, h=nTest)
liveness.ses.mape = mape(liveness.train.ses$mean, liveness.test)
plot(liveness.train.ses, main=paste("liveness", "SES"), sub=paste("MAPE:",round(liveness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#SARIMA
liveness.sarima.model <- auto.arima(liveness.train)
liveness.sarima <- forecast(liveness.sarima.model, h=nTest)
liveness.sarima.mape = mape(liveness.sarima$mean, liveness.test)
plot(liveness.sarima, main=paste("liveness", "SARIMA"), sub=paste("MAPE:", round(liveness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(liveness)

#Neural Network
liveness.train.nn <- elm(liveness.train)
liveness.nn.forecast <- forecast(liveness.train.nn, h=nTest)
liveness.nn.mape = mape(liveness.nn.forecast$mean, liveness.test)
plot(liveness, main=paste("liveness", "NN"), sub=paste("MAPE:", round(liveness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(liveness.nn.forecast$mean, start=1921+nTrain), col=2)

Loudness

loudness = ts(as.numeric(feature_quantiles[,"loudness"]) , start=1921)
loudness.train = subset(loudness, start=1, end=nTrain)
loudness.test = subset(loudness, start = (nTrain+1), end =(nTrain+nTest))

#SES
loudness.train.ses <- ses(loudness.train, h=nTest)
loudness.ses.mape = mape(loudness.train.ses$mean, loudness.test)
plot(loudness.train.ses, main=paste("loudness", "SES"), sub=paste("MAPE:",round(loudness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#SARIMA
loudness.sarima.model <- auto.arima(loudness.train)
loudness.sarima <- forecast(loudness.sarima.model, h=nTest)
loudness.sarima.mape = mape(loudness.sarima$mean, loudness.test)
plot(loudness.sarima, main=paste("loudness", "SARIMA"), sub=paste("MAPE:", round(loudness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(loudness)

#Neural Network
loudness.train.nn <- elm(loudness.train)
loudness.nn.forecast <- forecast(loudness.train.nn, h=nTest)
loudness.nn.mape = mape(loudness.nn.forecast$mean, loudness.test)
plot(loudness, main=paste("loudness", "NN"), sub=paste("MAPE:", round(loudness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(loudness.nn.forecast$mean, start=1921+nTrain), col=2)

Speechiness

speechiness = ts(as.numeric(feature_quantiles[,"speechiness"]) , start=1921)
speechiness.train = subset(speechiness, start=1, end=nTrain)
speechiness.test = subset(speechiness, start = (nTrain+1), end =(nTrain+nTest))

#SES
speechiness.train.ses <- ses(speechiness.train, h=nTest)
speechiness.ses.mape = mape(speechiness.train.ses$mean, speechiness.test)
plot(speechiness.train.ses, main=paste("speechiness", "SES"), sub=paste("MAPE:",round(speechiness.ses.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(speechiness)

#SARIMA
speechiness.sarima.model <- auto.arima(speechiness.train)
speechiness.sarima <- forecast(speechiness.sarima.model, h=nTest)
speechiness.sarima.mape = mape(speechiness.sarima$mean, speechiness.test)
plot(speechiness.sarima, main=paste("speechiness", "SARIMA"), sub=paste("MAPE:", round(speechiness.sarima.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(speechiness)

#Neural Network
speechiness.train.nn <- elm(speechiness.train)
speechiness.nn.forecast <- forecast(speechiness.train.nn, h=nTest)
speechiness.nn.mape = mape(speechiness.nn.forecast$mean, speechiness.test)
plot(speechiness, main=paste("speechiness", "NN"), sub=paste("MAPE:", round(speechiness.nn.mape*100, 3), "%"), xlab="Year", ylab="Median Value")
lines(ts(speechiness.nn.forecast$mean, start=1921+nTrain), col=2)

Creating linear models for predicting popular features based on their current feature

#Energy
best_energy = ts(as.numeric(best_feature_quantiles[,"energy"]), start=1921, end=2020) 
best_energy.train = subset(best_energy, start = 1, end = nTrain)
best_energy.test = subset(best_energy, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_energy.train, overall = energy.train)

best = best_energy.train
overall = energy.train

energy.lm <- lm(best ~ overall)
test_df <- data.frame(overall = energy.test)
best_energy.test.predict <- predict(energy.lm, newdata=test_df, interval="prediction")
best_energy_test_predict.mape <- mape(best_energy.test[1:length(best_energy.test)], best_energy.test.predict[,1])
plot(best_energy,xlab="Year", ylab="Energy", main=paste("Energy Population Prediction"), sub = paste("MAPE: ", round(best_energy_test_predict.mape*100, 3), "%"))
lines(ts(best_energy.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Duration
best_duration = ts(as.numeric(best_feature_quantiles[,"duration_ms"]), start=1921, end=2020) 
best_duration.train = subset(best_duration, start = 1, end = nTrain)
best_duration.test = subset(best_duration, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_duration.train, overall = duration.train)

best = best_duration.train
overall = duration.train

duration.lm <- lm(best ~ overall)
test_df <- data.frame(overall = duration.test)
best_duration.test.predict <- predict(duration.lm, newdata=test_df, interval="prediction")
best_duration_test.predict.mape <- mape(best_duration.test[1:length(best_duration.test)], best_duration.test.predict[,1])
plot(best_duration,xlab="Year", ylab="duration", main=paste("duration Popularity Prediction"), sub = paste("MAPE: ", round(best_duration_test.predict.mape*100, 3), "%"))
lines(ts(best_duration.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Accousticness
best_acousticness = ts(as.numeric(best_feature_quantiles[,"acousticness"]), start=1921, end=2020) 
best_acousticness.train = subset(best_acousticness, start = 1, end = nTrain)
best_acousticness.test = subset(best_acousticness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_acousticness.train, overall = acousticness.train)

best = best_acousticness.train
overall = acousticness.train

acousticness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = acousticness.test)
best_acousticness.test.predict <- predict(acousticness.lm, newdata=test_df, interval="prediction")
best_acousticness_test.predict.mape <- mape(best_acousticness.test[1:length(best_acousticness.test)], best_acousticness.test.predict[,1])
plot(best_acousticness,xlab="Year", ylab="acousticness", main=paste("acousticness Popularity Prediction"), sub = paste("MAPE: ", round(best_acousticness_test.predict.mape*100, 3), "%"))
lines(ts(best_acousticness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Valence
best_valence = ts(as.numeric(best_feature_quantiles[,"valence"]), start=1921, end=2020) 
best_valence.train = subset(best_valence, start = 1, end = nTrain)
best_valence.test = subset(best_valence, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_valence.train, overall = valence.train)

best = best_valence.train
overall = valence.train

valence.lm <- lm(best ~ overall)
test_df <- data.frame(overall = valence.test)
best_valence.test.predict <- predict(valence.lm, newdata=test_df, interval="prediction")
best_valence_test.predict.mape <- mape(best_valence.test[1:length(best_valence.test)], best_valence.test.predict[,1])
plot(best_valence,xlab="Year", ylab="Valence", main=paste("Valence Popularity Prediction"), sub = paste("MAPE: ", round(best_valence_test.predict.mape*100, 3), "%"))
lines(ts(best_valence.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Tempo
best_tempo = ts(as.numeric(best_feature_quantiles[,"tempo"]), start=1921, end=2020) 
best_tempo.train = subset(best_tempo, start = 1, end = nTrain)
best_tempo.test = subset(best_tempo, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_tempo.train, overall = tempo.train)

best = best_tempo.train
overall = tempo.train

tempo.lm <- lm(best ~ overall)
test_df <- data.frame(overall = tempo.test)
best_tempo.test.predict <- predict(tempo.lm, newdata=test_df, interval="prediction")
best_tempo_test.predict.mape <- mape(best_tempo.test[1:length(best_tempo.test)], best_tempo.test.predict[,1])
plot(best_tempo,xlab="Year", ylab="tempo", main=paste("tempo Popularity Prediction"), sub = paste("MAPE: ", round(best_tempo_test.predict.mape*100, 3), "%"))
lines(ts(best_tempo.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Liveness
best_liveness = ts(as.numeric(best_feature_quantiles[,"liveness"]), start=1921, end=2020) 
best_liveness.train = subset(best_liveness, start = 1, end = nTrain)
best_liveness.test = subset(best_liveness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_liveness.train, overall = liveness.train)

best = best_liveness.train
overall = liveness.train

liveness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = liveness.test)
best_liveness.test.predict <- predict(liveness.lm, newdata=test_df, interval="prediction")
best_liveness_test.predict.mape <- mape(best_liveness.test[1:length(best_liveness.test)], best_liveness.test.predict[,1])
plot(best_liveness,xlab="Year", ylab="liveness", main=paste("liveness Popularity Prediction"), sub = paste("MAPE: ", round(best_liveness_test.predict.mape*100, 3), "%"))
lines(ts(best_liveness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Loudness
best_loudness = ts(as.numeric(best_feature_quantiles[,"loudness"]), start=1921, end=2020) 
best_loudness.train = subset(best_loudness, start = 1, end = nTrain)
best_loudness.test = subset(best_loudness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_loudness.train, overall = loudness.train)

best = best_loudness.train
overall = loudness.train

loudness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = loudness.test)
best_loudness.test.predict <- predict(loudness.lm, newdata=test_df, interval="prediction")
best_loudness_test.predict.mape <- mape(best_loudness.test[1:length(best_loudness.test)], best_loudness.test.predict[,1])
plot(best_loudness,xlab="Year", ylab="loudness", main=paste("loudness Popularity Prediction"), sub = paste("MAPE: ", round(best_loudness_test.predict.mape*100, 3), "%"))
lines(ts(best_loudness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

#Speechiness
best_speechiness = ts(as.numeric(best_feature_quantiles[,"speechiness"]), start=1921, end=2020) 
best_speechiness.train = subset(best_speechiness, start = 1, end = nTrain)
best_speechiness.test = subset(best_speechiness, start = nTrain+1, end = nTrain+nTest)
train_df = data.frame(best = best_speechiness.train, overall = speechiness.train)

best = best_speechiness.train
overall = speechiness.train

speechiness.lm <- lm(best ~ overall)
test_df <- data.frame(overall = speechiness.test)
best_speechiness.test.predict <- predict(speechiness.lm, newdata=test_df, interval="prediction")
best_speechiness_test.predict.mape <- mape(best_speechiness.test.predict[,1], best_speechiness.test[1:length(best_speechiness.test)])
plot(best_speechiness,xlab="Year", ylab="speechiness", main=paste("speechiness Popularity Prediction"), sub = paste("MAPE: ", round(best_speechiness_test.predict.mape*100, 3), "%"))
lines(ts(best_speechiness.test.predict[,1], start = 1921+nTrain+1, end = 2020), col = 2)

Predict features for the overall market trends for the next five years. We will only use models that had the low enough MAPEs (10% chosen as a threshold) based on potential timeseries models.

Energy: NN Model

energy.nn <- elm(energy)
energy.nn.predict <- forecast(energy.nn, h=nTest)
plot(energy.nn.predict, main=paste("Energy", "NN"), xlab="Year", ylab="Median Value")
lines(energy)

Duration: NN Model

duration.nn <- elm(duration)
duration.nn.predict <- forecast(duration.nn, h=nTest)
plot(duration.nn.predict, main=paste("Duration", "NN"), xlab="Year", ylab="Median Value")

Tempo: SARIMA Model

tempo.sarima <- auto.arima(tempo)
tempo.sarima.predict <- forecast(tempo.sarima, h=nTest)
plot(tempo.sarima.predict, main=paste("Tempo", "SARIMA"), xlab="Year", ylab="Median Value")
lines(tempo)

Valence: SARIMA Model

valence.sarima <- auto.arima(valence)
valence.sarima.predict <- forecast(valence.sarima, h=nTest)
plot(valence.sarima.predict, main=paste("Valence", "SARIMA"), xlab="Year", ylab="Median Value")
lines(valence)

Liveness: SARIMA Model

liveness.sarima <- auto.arima(liveness)
liveness.sarima.predict <- forecast(liveness.sarima, h=nTest)
plot(liveness.sarima.predict, main=paste("liveness", "SARIMA"), xlab="Year", ylab="Median Value")
lines(liveness)

Loudness: SES Model

loudness.ses <- ses(loudness, h=nTest)
plot(loudness.ses, main=paste("loudness", "SES"), xlab="Year", ylab="Median Value")
lines(loudness)

Now predicting the median values for the most popular songs for this set of features.

Energy

best = best_energy
overall = energy


energy.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = energy.nn.predict$mean)
energy.predict <- predict(energy.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(energy, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Energy Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(energy.nn.predict$mean, 3), collapse=",") ))
lines(best_energy, col = 2)

lines(ts(energy.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Duration

best = best_duration
overall = duration

duration.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = duration.nn.predict$mean)
duration.predict <- predict(duration.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(duration, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Duration Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(duration.nn.predict$mean), collapse=",") ))
lines(best_duration, col = 2)

lines(ts(duration.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Tempo

best = best_tempo
overall = tempo

tempo.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = tempo.sarima.predict$mean)
tempo.predict <- predict(tempo.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(tempo, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's Tempo Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(tempo.sarima.predict$mean), collapse=",") ))
lines(best_tempo, col = 2)

lines(ts(tempo.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Valence

best = best_valence
overall = valence

valence.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = valence.sarima.predict$mean)
valence.predict <- predict(valence.lm, newdata=predict_df, interval="prediction")


smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(valence, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's valence Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(valence.sarima.predict$mean,3), collapse=",") ))
lines(best_valence, col = 2)

lines(ts(valence.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Liveness

best = best_liveness
overall = liveness

liveness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = liveness.sarima.predict$mean)
liveness.predict <- predict(liveness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(liveness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's liveness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(liveness.sarima.predict$mean,3), collapse=",") ))
lines(best_liveness, col = 2)

lines(ts(liveness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("topright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))

Loudness

best = best_loudness
overall = loudness

loudness.lm <- lm(best ~ overall)
predict_df <- data.frame(overall = loudness.ses$mean)
loudness.predict <- predict(loudness.lm, newdata=predict_df, interval="prediction")

smallest = min(min(best), min(overall))
biggest = max(max(best), max(overall))
plot(loudness, xlab="Year", ylab="Median Value", ylim=c(smallest, biggest), main="Most Popular Song's loudness Feature Prediction", sub=paste("Projected median feature value for most popular songs:", paste(round(loudness.ses$mean,1), collapse=",") ))
lines(best_loudness, col = 2)

lines(ts(loudness.predict[,1], start=2020+1, end=2020+nTest), col = 3)
legend("bottomright", lty=1, col=c(1,2,3), legend=c("Overall", "Most Popular", "Most Popular Projection"))